In [1]:
%load_ext autoreload
In [2]:
#load in modules
%autoreload 2
from wildfireassessment.ops import * #my package
import rasterio
from rasterio.transform import xy
from rasterio.plot import show, show_hist
from rasterio.mask import mask
import matplotlib.pyplot as plt
from pathlib import Path
import os
import earthpy as et
import earthpy.plot as ep
import numpy as np
import geopandas as gpd
from fiona.crs import from_epsg
from pyproj import Proj, transform
from functools import partial
import pyproj
from shapely.ops import transform
from rasterstats import zonal_stats
import itertools
%matplotlib inline
# from osgeo import gdal, osr, ogr, gdal_array
# from osgeo.gdalconst import *
Load WorldView images that are clipped to AOI¶
In [3]:
#read in filepaths for data
filepath_post = Path("./data/Paradise/post")
filepath_pre = Path("./data/Paradise/pre")
#WorldView Post/Pre
fps_wv_post = list(filepath_post.glob("2*_clip.tif"))
fps_wv_pre = list(filepath_pre.glob("2*_clip.tif"))
#Sent2 Post/Pre
fp_sent2_post = filepath_post / "B08_post_clipped_byte_scaled.tif"
fp_sent2_pre = filepath_pre / "B08_pre_clipped_byte_scaled.tif"
In [8]:
raster_src_post, rgb_post = readRGBImg(fps_wv_post[5])
raster_src_pre, rgb_pre = readRGBImg(fps_wv_pre[5])
In [9]:
raster_src_post_b08, b08_post = readOneImg(fp_sent2_post)
raster_src_pre_b08, b08_pre = readOneImg(fp_sent2_pre)
In [7]:
!gdalinfo data/Paradise/post/B08_post_clipped_byte_scaled.tif
Clipped images¶
In [10]:
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(15, 10))
im = ax1.imshow(rgb_pre)
ax1.set_title("Pre-fire")
im = ax2.imshow(rgb_post)
ax2.set_title("Post-fire")
plt.show()
In [11]:
x_r_pre = rgb_pre[:,:,0].ravel().reshape(-1, 1)
x_g_pre = rgb_pre[:,:,1].ravel().reshape(-1, 1)
x_b_pre = rgb_pre[:,:,2].ravel().reshape(-1, 1)
x_pre = np.hstack((x_r_pre, x_g_pre, x_b_pre))
x_r_pre, x_g_pre, x_b_pre = None, None, None
y_r_post = rgb_post[:,:,0].ravel().reshape(-1, 1)
y_g_post = rgb_post[:,:,1].ravel().reshape(-1, 1)
y_b_post = rgb_post[:,:,2].ravel().reshape(-1, 1)
y_post = np.hstack((y_r_post, y_g_post, y_b_post))
y_r_post, y_g_post, y_b_post = None, None, None
In [12]:
n_bins = 20
colors = ['red', 'green', 'blue']
fig, axs = plt.subplots(1, 2, sharey=True, tight_layout=True, figsize=(15, 10))
axs[0].hist(x_pre, n_bins, alpha=0.5, color=colors, label=colors)
axs[1].hist(y_post, n_bins, alpha=0.5, color=colors, label=colors)
axs[0].legend(prop={'size': 10})
axs[1].legend(prop={'size': 10})
axs[0].set_title("Pre-fire")
axs[1].set_title("Post-fire")
plt.show()
MinMax Histogram Stretch on Sent2 images¶
MinMax histogram stretch on values on Sent2 images. The maximum was determined in properties of pre-fire image. This is done since DigitalGlobe images were adjusted according to their Dynamic Range Adjustment (no specifics on how this histogram stretch is done).
In [13]:
#command in gdal translate
#!gdal_translate -scale 0 7210 0 255 -ot Byte data/Paradise/post/B08_post_clipped.tif data/Paradise/post/B08_post_clipped_byte_scaled.tif
In [14]:
#command in gdal translate
#!gdal_translate -scale 0 7210 0 255 -ot Byte data/Paradise/pre/B08_pre_clipped.tif data/Paradise/pre/B08_pre_clipped_byte_scaled.tif
In [15]:
bbox_post = box(raster_src_post.bounds.left, raster_src_post.bounds.bottom, raster_src_post.bounds.right, raster_src_post.bounds.top)
bbox_pre = box(raster_src_pre.bounds.left, raster_src_pre.bounds.bottom, raster_src_pre.bounds.right, raster_src_pre.bounds.top)
In [16]:
out_img_post_b08, out_img_transform_post_b08, out_meta_post = clipImg(bbox_post, raster_src_post_b08, proj=4326)
out_img_pre_b08, out_img_transform_pre_b08 , out_meta_pre = clipImg(bbox_pre, raster_src_pre_b08, proj=4326)
In [17]:
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(15, 10))
im = ax1.imshow(out_img_pre_b08[0], cmap='gray')
ax1.set_title("Pre-fire")
im = ax2.imshow(out_img_post_b08[0], cmap='gray')
ax2.set_title("Post-fire")
plt.show()
In [18]:
x_nir_pre = out_img_pre_b08[0].ravel()
y_nir_post = out_img_post_b08[0].ravel()
In [19]:
n_bins = 20
colors = ['red']
fig, axs = plt.subplots(1, 2, sharey=True, tight_layout=True, figsize=(8, 5))
axs[0].hist(x_nir_pre, n_bins, alpha=0.5, color=colors, label='nir')
axs[1].hist(y_nir_post, n_bins, alpha=0.5, color=colors, label='nir')
axs[0].legend(prop={'size': 10})
axs[1].legend(prop={'size': 10})
axs[0].set_title("Pre-fire")
axs[1].set_title("Post-fire")
plt.show()
Segment image¶
In [24]:
indices = chunkImageIndices(rgb_post)
indices
Out[24]:
In [26]:
coord_ind = (9792-1024, 9792+1024, 9792-1024, 9792+1024)
In [59]:
from skimage.segmentation import slic, felzenszwalb
from skimage.segmentation import mark_boundaries
from skimage.measure import label, regionprops
segments_slic = slic(rgb_pre[coord_ind[0]:coord_ind[1], coord_ind[2]:coord_ind[3], :], n_segments=1000, compactness=15)
#segments_slic = slic(rgb_pre[indices[1][0]:indices[1][1], indices[1][2]:indices[1][3], :], n_segments=5000, compactness=10)
#segments_slic = np.add(np.ones(segments_slic.shape), segments_slic).astype(np.uint16)
print('SLIC number of segments: {}'.format(len(np.unique(segments_slic))))
#segments_fz = felzenszwalb(out_img_pre[0:chunksize[0], 0:chunksize[1], :], scale=10, sigma=0.5, min_size=20)
#print('Felz number of segments: {}'.format(len(np.unique(segments_fz))))
In [60]:
boundary_rgb_post = mark_boundaries(rgb_post[coord_ind[0]:coord_ind[1], coord_ind[2]:coord_ind[3], :], segments_slic, outline_color=(1, 1, 0), mode='thick')
bouncary_rgb_pre = mark_boundaries(rgb_pre[coord_ind[0]:coord_ind[1], coord_ind[2]:coord_ind[3], :], segments_slic, outline_color=(1, 1, 0), mode='thick')
# boundary_rgb_post = mark_boundaries(rgb_post[indices[1][0]:indices[1][1], indices[1][2]:indices[1][3], :], segments_slic, outline_color=(1, 1, 0), mode='thick')
# bouncary_rgb_pre = mark_boundaries(rgb_pre[indices[1][0]:indices[1][1], indices[1][2]:indices[1][3], :], segments_slic, outline_color=(1, 1, 0), mode='thick')
#boundary_rgb_post = mark_boundaries(out_img_pre[0:chunksize[0], 0:chunksize[1], :], segments_fz, outline_color=(1, 1, 0), mode='thick')
fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(15, 30))
im = ax1.imshow(bouncary_rgb_pre)
ax1.set_title("Pre-fire")
im = ax2.imshow(boundary_rgb_post)
ax2.set_title("Post-fire")
plt.show()
Vectorize segments with spectral vals¶
In [29]:
from affine import Affine
In [61]:
(indices[0][1]*indices[0][3]//(2048*2048))*1000
Out[61]:
In [42]:
indices[1]
Out[42]:
In [37]:
raster_src_pre.transform
Out[37]:
In [44]:
raster_src_pre.index(raster_src_pre.transform[2], raster_src_pre.transform[5])
Out[44]:
In [41]:
Affine.translation(1.0, 1.0) * raster_src_pre.transform
Out[41]:
In [24]:
raster_src_pre.height
Out[24]:
In [50]:
tup = raster_src_pre.xy(8168, 0)
print(tup)
tup_origin = raster_src_pre.xy(0,0)
print(tup_origin)
In [51]:
abs(tup[0] - tup_origin[0])
Out[51]:
In [52]:
abs(tup[1] - tup_origin[1])
Out[52]:
In [57]:
x_trans = abs(tup[0] - tup_origin[0])
y_trans = abs(tup[1] - tup_origin[1])
Affine.translation(x_trans, -y_trans) * raster_src_pre.transform
Out[57]:
In [22]:
raster_src_post.crs
Out[22]:
In [58]:
import shapely
from shapely.wkt import loads
from shapely.geometry import Polygon, mapping, shape, MultiPolygon, box
from rasterio import features
import geopandas as gpd
import pandas as pd
shapes_list = [{'seg_index': int(v), 'geometry': loads(shape(g).wkt)} for g, v in features.shapes(segments_slic.astype(np.uint16), mask=None, transform=Affine.translation(x_trans, -y_trans) * raster_src_pre.transform)]
gdf = gpd.GeoDataFrame(shapes_list)
gdf.crs = {'init': 'EPSG:4326'}
gdf.insert(2, 'crs', str(raster_src_post.crs))
gdf.head()
Out[58]:
In [59]:
gdf.to_file("segments.shp")
In [41]:
## use regionprops onn segments to extract properties for dataframe
#post
regions_red = regionprops(segments_slic, rgb_post[0:chunksize[0], 0:chunksize[1],0])
regions_blue = regionprops(segments_slic, rgb_post[0:chunksize[0], 0:chunksize[1],1])
regions_green = regionprops(segments_slic, rgb_post[0:chunksize[0], 0:chunksize[1],2])
#pre
regions_red_pre = regionprops(segments_slic, rgb_pre[0:chunksize[0], 0:chunksize[1],0])
regions_blue_pre = regionprops(segments_slic, rgb_pre[0:chunksize[0], 0:chunksize[1],1])
regions_green_pre = regionprops(segments_slic, rgb_pre[0:chunksize[0], 0:chunksize[1],2])
In [42]:
gdf = gdf[gdf['seg_index'] != 0]
In [43]:
project = partial(
pyproj.transform,
pyproj.Proj(init='epsg:4326'), # source coordinate system
pyproj.Proj(init='epsg:32610')) # destination coordinate system
# g1 = gdf.loc[gdf['seg_index']==1000, 'geometry'].iloc[0]
# g2 = transform(project, g1) # apply projection
In [44]:
# create a dataframe with spectral values as column names , index value "seg_index" for merging geodataframe with Polygons
region_spectrals = []
for i in range(len(regions_red)):
seg_label = regions_red[i].label
g1 = gdf.loc[gdf['seg_index']==seg_label, 'geometry'].iloc[0]
g2 = transform(project, g1)
nir_zone_post = zonal_stats(g2, out_img_post_b08[0], affine=out_img_transform_post_b08, stats='mean', nodata=-999)
nir_zone_pre = zonal_stats(g2, out_img_pre_b08[0], affine=out_img_transform_pre_b08, stats='mean', nodata=-999)
dict_seg = {'seg_index': regions_red[i].label,
'red_value': regions_red[i].mean_intensity,
'blue_value': regions_blue[i].mean_intensity,
'green_value': regions_green[i].mean_intensity,
'nir_value': nir_zone_post[0]['mean'],
'red_value_pre': regions_red_pre[i].mean_intensity,
'blue_value_pre': regions_blue_pre[i].mean_intensity,
'green_value_pre': regions_green_pre[i].mean_intensity,
'nir_value_pre': nir_zone_pre[0]['mean'],
'area_m': regions_red[i].area * 0.31} #area in meters
region_spectrals.append(dict_seg)
In [46]:
df = pd.DataFrame(region_spectrals)
df.head()
Out[46]:
In [47]:
# merge on seg_index
gdf2 = pd.merge(gdf, df, on='seg_index')
# #add in these two columns for labelling later
# gdf2['land_class'] = np.nan
# gdf2['burn_class'] = np.nan
gdf2.head()
Out[47]:
Compute SIs¶
$$SI= \frac{Band_1-Band_2}{Band_1+Band_2}$$$$\Delta SI= SI_{pre}-SI_{post}$$In [48]:
def computeSI(b1, b2):
return (b1-b2)/(b1+b2)
def changedSI(SI_pre, SI_post):
return SI_pre - SI_post
In [49]:
# convert the keys to list
post_keys = ['blue_value', 'green_value', 'red_value', 'nir_value']
pre_keys = ['blue_value_pre', 'green_value_pre', 'red_value_pre', 'nir_value_pre']
In [50]:
perm = itertools.permutations(post_keys, 2)
perm_pre = itertools.permutations(pre_keys, 2)
In [51]:
post_SI_keys = []
for i, tup in enumerate(perm):
b1 = tup[0]
b2 = tup[1]
col_name = "SI_" + b1[0] + b2[0] + "_post"
post_SI_keys.append(col_name)
print(col_name)
gdf2[col_name] = computeSI(gdf2[b1], gdf2[b2])
gdf2.head()
Out[51]:
In [52]:
pre_SI_keys = []
for i, tup in enumerate(perm_pre):
b1 = tup[0]
b2 = tup[1]
col_name = "SI_" + b1[0] + b2[0] + "_pre"
pre_SI_keys.append(col_name)
print(col_name)
gdf2[col_name] = computeSI(gdf2[b1], gdf2[b2])
gdf2.head()
Out[52]:
compute dSIs
In [53]:
SI_keys = pre_SI_keys + post_SI_keys
SI_keys
Out[53]:
In [54]:
band_combos = list(set([key.split('_')[1] for key in post_SI_keys]))
SI_combos = [ tuple([s for s in SI_keys if bcombostr in s]) for bcombostr in band_combos]
SI_combos
Out[54]:
In [55]:
# now, add dSIs to dataframe
for i, tup in enumerate(SI_combos):
SI_post = tup[0]
SI_pre = tup[1]
col_name = "dSI_" + SI_post.split('_')[1]
gdf2[col_name] = changedSI(gdf2[SI_pre], gdf2[SI_post])
gdf2.head()
Out[55]:
In [56]:
gdf2['land_class'] = np.nan
gdf2['burn_class'] = np.nan
gdf2.head()
Out[56]:
In [57]:
gdf2.to_file("./data/segments.geojson", driver='GeoJSON')
write data to files¶
In [7]:
fps_wv_post[5:6]
Out[7]:
In [ ]:
writeDatasets(fps_wv_post[5:6], fps_wv_post[5:6], fp_sent2_post, fp_sent2_pre)
#writeDatasets(fps_wv_post[8:], fps_wv_pre[8:], fp_sent2_post, fp_sent2_pre)
In [ ]: